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The dependencies on strain and oxygen vacancies of the ferroelectric polarization and the weak 
ferromagnetic magnetization in the multiferroic material bismuth ferrite, BiFeOs, are investigated 
using first principles density functional theory calculations. The electric polarization is found to 
be rather independent of strain, in striking contrast to most conventional perovskite ferroelectrics. 
It is also not significantly affected by oxygen vacancies, or by the combined presence of strain 
and oxygen vacancies. The magnetization is also unaffected by strain, however the incorporation of 
oxygen vacancies can alter the magnetization slightly, and also leads to the formation of Fe 2+ . These 
results are discussed in light of recent experiments on epitaxial films of BiFeC>3 which reported a 
strong thickness dependence of both magnetization and polarization. 

PACS numbers: 71.15.Mb, 75.70.Ak, 77.80.-e 
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I. INTRODUCTION 

Materials that simultaneously show electric and mag- 
netic order are currently gaining more and more atten- 
tion. This is partly due to the fact that such multiferroics 
are promising materials for new types of multifunctional 
device applications, but also because of the interesting 
physics found in this class of materials. For example a 
strong coupling between ferroelectric and antiferromag- 
netic domain walls has been found in YMnOsji in or- 
thorhombic TbMnC>3 and TbMn204 the ferroelectric po- 
larization can be reoriented by a magnetic fieldj 2 ^ and 
ferromagnetic ordering can be "switched on" by an elec- 
tric field in hexagonal HoMnOs^ 

Although magnetoelectric materials have been known 
for a long time^ recent progress in thin-film growth and 
other sample preparation techniques contributed consid- 
erably to their renaissance. By using techniques such as 
pulsed laser deposition (PLD), chemical vapor deposition 
(CVD), or molecular beam epitaxy (MBE), many mate- 
rials can nowadays be prepared as high quality epitaxial 
thin films. One advantage of these techniques is the pos- 
sibility to stabilize otherwise metastable structures or to 
tune material properties by varying the lattice mismatch 
between the film and the substrate, thereby introducing 
epitaxial strain in the thin film material. 

Indeed, the effect of strain on the ferroelectric prop- 
erties of conventional ferroelectric materials is also a 
topic of current interest. Strain effects can lead to a 
substantial increase of the spontaneous polarization and 
Curie temperature^ and even drive paraelectric materi- 
als (such as SrTiC-3) into the ferroelectric phased Since 
the mechanism for ferroelectricity in multiferroic materi- 
als is often different from that in conventional perovskite 
ferroelectrics r& the question arises of whether similar 
strain effects will be observed in multiferroics. Magnetic 
properties can also be strongly affected by strain, mainly 
due to large changes in anisotropy. 10 Strain can also affect 
the saturation magnetization and Curie temperature. 11 



First principles density functional theory (DFT) cal- 
culations (see e.g. Ref. ll2l) play a crucial role in studying 
the influence of strain on ferroelectric propertiesiiii^ii^ 
Epitaxial strain can be introduced straightforwardly in 
DFT studies by fixing the lattice parameters in the di- 
rections corresponding to the lateral dimension of the 
substrate and allowing the system to relax in the per- 
pendicular direction. This makes it possible to clearly 
distinguish between the effect of strain and other influ- 
ences present in real thin-film samples, such as interface 
effects or various types of defects. Such information can 
then be used to optimize the properties of the thin-film 
material. 

In this work we study the influence of strain on the 
electric polarization and magnetization of multiferroic 
bismuth ferrite, BiFe03. Bismuth ferrite crystallizes 
in a rhombohedrally distorted perovskite structure with 
space group i?3cji& where all ions are displaced along 
the [111] direction relative to the ideal centrosymmet- 
ric positions, and the oxygen octahedra surrounding the 
transition metal cations are rotated around this axis, 
alternately clockwise and counterclockwise. The mag- 
netic order is essentially G-type antiferromagnetioii but 
in addition, the direction of the magnetic moments in the 
bulk rotates with a long wavelength of 620 A(Ref. GJ). 
We have recently shown that if this spiral spin struc- 
ture is suppressed, the system shows weak ferromag- 
netismMi with the magnetic moments oriented perpen- 
dicular to the rhombohedral axis and a slight canting of 
these magnetic moments resulting in a small macroscopic 
magnetization^ 

Ferroelectric hysteresis loops have been measured for 
BiFeOs but the experimental determination of an exact 
value for the spontaneous polarization P s is difficult due 
to large leakage currents*^ Several values are reported 
in the literature, summarized in Table I of Ref. 1 2 J . A 
first principles calculation of the spontaneous polariza- 
tion of bulk BiFeOs results in a value of P s ~ 95 /iC/cm 2 
(Ref. l22r> : recent measurements for epitaxial films grown 
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on SrTiC>3 agree well with this value ^^S^ These exper- 
iments also show a strong dependence of both magnetiza- 
tion and ferroelectric polarization on the film thickness. 23 
One likely explanation for this thickness dependence is 
the increase in strain with decreasing film thickness; an- 
other is a change in the concentration of defects such as 
oxygen vacancies. In this work we systematically exam- 
ine both possibilities using first principles DFT calcula- 
tions. 

Our main result is that the strain dependence of the 
ferroelectric polarization in BiFe03 is rather weak com- 
pared with conventional ferroelectric materials and that 
it cannot explain the variation of the polarization re- 
ported for the thin film samples. The same is true for the 
magnetization which also shows only weak strain depen- 
dence. The weak strain dependence of the polarization is 
due to a very stable ionic configuration in BiFe03, which 
manifests itself in only small changes of the relative ionic 
positions when the lattice is strained. 

The high stability of the ferroelectric configuration also 
leads to a negligible dependence of the electric polariza- 
tion on the oxygen vacancy concentration. In contrast, 
we find that the magnetization of BiFe03 is affected by 
the presence of oxygen vacancies but the changes are not 
very systematic and depend on the precise position of 
the oxygen vacancy. The presence of oxygen vacancies 
in all cases leads to the formation of Fe 2+ ions, which 
can be identified unequivocally in the partial densities of 
states although the actual charge differences between the 
different Fe sites are very small. 



II. METHOD 

For this work we use first principles density func- 
tional theory (see e.g. Ref. ll2T) within the projector aug- 
mented wave (PAW) method— as implemented in the 
Vienna Ab-initio Simulation Package (VASP)^*2&. We 
include 15 valence electrons for Bi (5<i 10 6s 2 6p 3 ), 14 for 
Fe (3p 6 3d 6 4s 2 ), and 6 for each oxygen (2s 2 2p 4 ), use an 
energy cutoff between 450-500 eV for the plane wave ex- 
pansion of the PAWs, a 4 x 4 x 4 Monkhorst Pack grid 
of k-pointS)S£ and the tetrahedron method with Blochl 
corrections for the Brillouin zone integrations^ For all 
structures (unless otherwise noted) we relax the ionic po- 
sitions while keeping the lattice vectors fixed until the 
Hellman-Feynman forces are less than 10~ 2 eV/A. For 
the calculation of the local densities of states at the Fe 
sites we use a sphere radius of 1.4 A. These values have 
been found to give good convergence of all quantities un- 
der consideration. 

To treat exchange and correlation effects we use both 
the local spin density approximation (LSDA) 1 - and the 
semiempirical LSDA+U method 3 ^, for a better descrip- 
tion of the localized transition metal d electrons. We 
have recently shown that using the LSDA+U method 
and a moderate value of U = 3 eV (and J = 1 eV) leads 
to a good description of the structural parameters and 



the ferroelectric polarization in BiFeOs^ Larger U val- 
ues shift the d bands further down in energy relative to 
the oxygen p states but have only a small effect on the 
structural and ferroelectric properties. U — 3 eV can be 
regarded as a lower limit of what is required to ensure 
the insulating character of BiFe03, and here we exclu- 
sively use this value (and J = 1 eV) in our LSDA+U 
calculations. We do not claim that these values necessar- 
ily would also lead to a good description of spectroscopic 
quantities such as e.g. photoemission spectra. 

There are two different LSDA+U approaches imple- 
mented in the VASP code, (i) the traditional LSDA+U 
approach of Anisimov, Liechtenstein and coworkers (in 
the so-called "fully localized limit" )£i, and (ii) a simpli- 
fied approach of Dudarev et alm^ where only the differ- 
ence U c ff = U — J enters. As shown in Appendix 1X1 the 
latter approach (ii) is identical to approach (i) when J 
— 0, so that the difference between these two approaches 
can be discussed in terms of a J-dependence. As pointed 
out above, the structural and ferroelectric properties do 
not depend strongly on the precise values of U and J and 
are therefore basically identical for both approaches. In 
contrast, the magnetization of BiFe03, which is due to 
a small canting of the mainly antiferromagnetically cou- 
pled magnetic moments of the Fe cationsjSI is strongly 
J-dependent. Although the absolute value of this canting 
(and therefore the macroscopic magnetization) depends 
strongly on J, the effects of strain and oxygen vacan- 
cies are basically independent of the actual value of J. 
Since the focus of the present paper is on these effects, 
we postpone the detailed analysis of the J dependence 
of the magnetization to a future publication and always 
present two data sets for the magnetization, one obtained 
using U = 3 eV and J = 1 eV, the second obtained us- 
ing U = 2 eV and J = eV (or equivalently approach 
(ii) with U C H = 2 eV). Our conclusions regarding the 
effects of strain and oxygen vacancies on the magnetiza- 
tion in BiFeOs apply independently to both data sets. 
We point out that a significant J dependence is only ob- 
served for the canting of the local magnetic moments of 
the Fe cations. The absolute values of these magnetic 
moments are rather independent of J, as long as J is 
varied within reasonable limits (J < 1.5 eV). 

For the calculation of the ferroelectric polarization we 
use the Berry-phase approach developed by Vanderbilt 
and King-SmithiS* 3 ^ 3 ^ In this theory the polarization 
of a periodic solid is represented by a three-dimensional 
lattice, and experimentally accessible polarization differ- 
ences are obtained by connecting two points of the "po- 
larization lattices" of the initial and final states, which 
can be transformed into each other through a continuous 
"path" of insulating intermediate states. In the present 
case special care must be taken in determining which 
points of the initial and final state polarization lattices 
have to be connected. First, because the value of the 
polarization difference in BiFeOs is comparable to the 
distance between neighboring points on the polarization 
lattices^ and in addition, because the direction of the 
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spontaneous polarization in the monoclinically strained 
structures and in the supercells containing oxygen vacan- 
cies is not symmetry restricted. In such cases the non- 
centrosymmetric distortions of the ionic positions in the 
corresponding systems are gradually reduced until an un- 
ambiguous connection can be made, i.e. the polarization 
of some intermediate states is explicitly calculated. 



III. RESULTS AND DISCUSSION 

A. Strain dependence of the electric polarization 
for (111) oriented films 

We first investigate the effect of strain corresponding 
to a (111) orientation of the substrate. This geome- 
try preserves the rhombohedral symmetry found in the 
unstrained system and the spontaneous polarization re- 
mains oriented along the [111] direction. We consider 
compressive strain as well as tensile strain. We fix the 
nearest neighbor distance a^ex between identical cations 
within the (111) plane to the values shown in Tableland 
then vary the out-of-plane lattice parameter Chex while 
relaxing all ionic positions^ The spontaneous electric 
polarization P s and magnetization per Fe cation M s are 
calculated for the relaxed value of the out-of-plane pa- 
rameter, c^ x , corresponding to the minimum of the total 
energy with fixed in-plane parameter ahox (and relaxed 
ionic coordinates). The results are summarized in Ta- 
ble The strain is defined as e = ah C x/&hcx.o — 1 where 
Q-hcx.o is the corresponding lattice constant for the un- 
strained system. For comparison we mention that the 
lattice constant of SrTi03, which is a commonly used 
substrate material, would lead to a compressive strain of 
e « -2 %. 

One can see that compressive epitaxial strain in the 
(111) plane leads to a reduction of the unit cell volume 
whereas tensile strain leads to a volume increase, i.e. the 
system does not behave like an ideal elastic medium. 
This result is similar to that found in Refs. |3y and [13| 
for BaTiC>3. The spontaneous polarization P s of BiFeC>3 
increases slightly with increasing compressive strain (4 % 
increase of the polarization for e = —3 %) but the ef- 
fect is much smaller than in the conventional ferroelec- 
tric systems BaTiC>3 and PbTiC>3. In BaTiC>3 a strain 
of e = — 1 % leads to an increase in the polarization of 
- 35 % (Ref.© whereas a similar strain in PbTiCK in- 
creases the polarization by 15-20 The largest strain 
effect on the polarization is seen in SrTiC"3 which is not 
ferroelectric in the unstrained state but develops a spon- 
taneous polarization in epitaxially strained films. 7 

We point out that the observed small change in spon- 
taneous polarization is consistent with the small changes 
in ionic displacements that result from the applied strain. 
If we calculate the expected change in polarization based 
on a simple point charge model as -P s (e) = P s (e = 
0) + y^y Ei Zi{Ri{e) -Ri(e= 0)) (here Z, is the charge 
associated with ion i (see below), Ri(e) is the correspond- 




FIG. 1: Polarization P s as function of strain e calculated 
within the Berry phase approach (full circles) and by a simple 
point charge model (see text) using the formal charges (open 
triangles) and Born-effective charges of the unstrained system 
(open squares). 



ing strain-dependent ionic position, V is the unit cell vol- 
ume (also strain dependent), and the sum extends over 
all ions in the unit cell) then the resulting change in po- 
larization for reasonable values of Zi is comparable to 
the change found by the direct calculation of P s using 
the Berry-phase approach. Fig. ^ shows the change in 
polarization as a function of strain, calculated with the 
Berry-phase approach as well as with the above formula 
using for Zi (i) the formal charges and (ii) the Born effec- 
tive charges (BECs) of the unstrained structured The 
weaker strain dependence of the spontaneous polarization 
in BiFeC>3 compared to other ferroelectrics can therefore 
be traced back to the weaker strain dependence of the 
ionic displacements in this material. 

Table [I] also shows the change in ionic displacements 
(compared to an ideal centrosymmetric reference struc- 
ture with the same lattice parameters) as a function of 
strain. The positions of the Bi cations are used as ref- 
erence and therefore the displacements of these ions are 
zero by definition. For a compressive strain of e = —3 % 
the displacements of the Fe cations change only by ~ 
17 % (LSDA) whereas the displacements of the oxygen 
atoms are nearly strain-independent. For comparison, an 
epitaxial strain of -2.28 % in the (001) plane of tetrago- 
nal BaTiC>3 leads to a change of the displacements of the 
Ti cations along [001] of ~ 52 %, and the displacements 
of the O anions in this case change by even more than 
100 %M 

There are three likely explanations for the weak strain 
dependence of the ionic displacements in BiFeC>3. The 
first is the general high stability of ferroelectrics with 
high Curie temperatures, large ionic displacements, and 
large energy differences between the ground state and 
the centrosymmetric reference structure. BiFe03 has a 
very high Curie temperature (Tc — 1123 K) and the en- 
ergy gain in the ground state R3c structure is ~ 0.25 eV 
per formula unit compared to the centrosymmetric R3c 
structure and ~ 1 eV compared to the cubic perovskite 
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TABLE I: Strain e, in-plane lattice parameter ahcx, relaxed out-of-plane lattice parameter c!^J x , volume V of the rhombo- 
hedral unit cell (containing two formula units), absolute and relative displacements of Fe and O ions along the [111] direc- 
tion (compared to an ideal centrosymmetric reference structure with the same lattice parameters, u,:(e) = Ri(e) — i?i,o(e), 
Aui(e) = («j(e)/uj(0)) — 1), spontaneous polarization P s , and magnetization M s for BiFe03 strained within the (111) plane. 
Upper values are obtained using the LSDA, lower values are obtained using the LSDA+U method. For the LSDA+U mag- 
netization, the first value refers to U = 2 eV, J = eV and the second value refers to U = 3 eV, J = 1 eV (see Sec. HTl for 
details) . 



e[%] 


ahcx [A] 




v [A 3 ] 


UFo [A] 


Am Pc [%] 


u [A] 


Au Q [%] 


P s [^C/cm 2 ] 


M e [mb/Fc] 


-3 


5.326 


13.68 


112.01 


0.294 


16.8 


0.527 


2.5 


102.8 


0.04 


-1 


5.431 


13.45 


114.53 


0.268 


6.2 


0.520 


1.2 


100.1 


0.05 





5.485 


13.31 


115.62 


0.252 





0.514 





98.9 


0.05 


+1 


5.541 


13.22 


117.16 


0.246 


-2.5 


0.517 


0.6 


98.4 


0.05 


+3 


5.655 


13.00 


120.01 


0.227 


-9.7 


0.520 


1.3 


97.9 


0.05 


-3 


5.343 


13.92 


114.71 


0.357 


21.8 


0.573 


6.2 


97.7 


0.02/0.11 





5.503 


13.48 


117.86 


0.293 





0.539 





94.0 


0.03/0.10 



structure). This explanation is consistent with the fact 
that the strain dependence of the polarization is already 
weaker in PbTi0 3 (Ref.Q T c = 763 K) than in BaTi0 3 
(Ref. [H T c = 400 K). If this explanation is valid, a 
similar strain independence of the electric polarization 
should be observable for LiNbC>3 which is isostructural 
to BiFeOa and has a even higher Curie temperature of 
1480 K (all Curie temperatures are taken from Ref. W^i . 
The second possible explanation is that the mechanism 
driving the ferroelectric distortion in BiFe03, namely the 
stereochemically active lone pairj2& is relatively inert to 
the changes in the lattice vectors caused by epitaxial 
strain. In conventional ferroelectric perovskites the ferro- 
electric distortion is stabilized by charge transfer from the 
oxygen into the unoccupied transition metal d orbitals,-- 
This charge transfer mechanism is probably more sen- 
sitive to small changes in bond lengths than the stereo- 
chemical activity of the Bi lone electron pair. In this case 
the study of the strain dependence of the electric polar- 
ization in the multiferroic BiMn03 would be of interest. 
The ferroelectricity in BiMnOs is driven by the Bi lone 
pair— but its "general stability" is not large (reported 
T c ~ 760 K, Ref. l4fj|) Thus a small (large) strain depen- 
dence in LiNbOs and a large (small) strain dependence 
in BiMnOs would confirm the "general stability" (lone 
pair) origin of the polarization stability. A third possible 
explanation for the weak strain dependence in BiFe03 is 
the special geometry of the oxygen octahedra in the Ric 
structure. In ferroelectrics like BaTiOs and PbTiC>3 all 
ions are displaced only along the polar direction whereas 
in BiFe03 (and also in LiNbOs) the oxygen octahedra 
are also rotated around this axis. The resulting geome- 
try of the oxygen cages surrounding the transition metal 
cations could be less favorable for an additional strain- 
induced displacement of the ions compared to the simpler 
geometry found in BaTiOs and PbTiC>3. Future studies 
will shed more light on this issue. 



The results obtained using the LSDA+U method are 
very similar to the LSDA results. Although the relative 
change in displacements seems to be slightly larger than 
within the LSDA (see Table P) , the relative changes in 
polarization are exactly the same in both cases. This 
reflects the fact, which has been already pointed out in 
Ref. I22I that the explicit treatment of electronic correla- 
tions within the LSDA+U method has only minor effect 
on the structural properties of this system. Apart from 
resulting in a slightly larger equilibrium volume, in better 
agreement with experimental data, the main effect is an 
improved description of the electronic structure resulting 
in a larger band gap and a stable insulating phase. 



B. Strain dependence of the polarization for (001) 
oriented films 

We now discuss the case of BiFeOs deposited on a (001) 
surface. In this case the epitaxial constraint of the cubic 
substrate enforces 90° angles of the lattice vectors within 
the (001) planes, in conflict with the rhombohedral dis- 
tortion of bulk BiFeOs . One can expect that the compe- 
tition between these two effects leads to a base-centered 
monoclinic structure and indeed monoclinic symmetry 
is found experimentally in epitaxial BiFeOs films de- 
posited on a (001) surface of SrTi03— The low sym- 
metry of the monoclinic structure renders a systematic 
computational investigation unfeasible for this type of 
epitaxial strain. We have therefore performed calcula- 
tions for two structures, one representative of very thin 
films (~ 100-200 nm), the other representative of thicker 
films (~ 400 nm), based on experimentally determined 
lattice parameters^ To account for the fact that the 
theoretically determined lattice parameters are usually 
slightly different from the experimental lattice parame- 
ters we have scaled all values accordingly, so that the 
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TABLE II: Experimental lattice parameters (in A) found in 
representative BiFeOs films— (line "Exp."), together with 
the corresponding lattice parameter a r h of rhombohedral bulk 
BiFeOs (Ref.©, and the values used in the calculations for 
the monoclinically strained films (line "Theo."). a and b cor- 
respond to the two in-plane directions and c to the out-of- 
plane direction of the monoclinic lattice. The values in the 
given form correspond to the lengths of the cube edges of the 
distorted cubic perovskite structure. Strain values e are also 
given. 





a r h/V2 


"thin films" 


"thick films" 


a/V2 b/V2 c 


a/s/2 b/y/2 c 


Exp. 
Theo. 


3.98 
3.89 


3.92 3.92 4.06 
3.83 3.83 3.97 


3.91 3.97 4.00 
3.82 3.88 3.91 


e 


-1.5% -1.5% +2% 


-2% -0.2% +0.4% 



TABLE III: Calculated absolute values of the polarization 
\P S \ and out-of-plane component (P 3 )[ooi] ( m /iC/cm 2 ) for 
the monoclinic structures. P CX p is the value measured in the 
[001] direction, 23 and should be compared to the calculated 
(P s )[ooi]- For the magnetization M 3 , the first value refers to 
U = 2 eV, J — eV and the second value refers to U = 3 eV, 
J — 1 eV (see Sec. [n] for details). 





\Ps\ 


(P>)[001] 


Pexp (Ref. 23J 


M B [/iB /Fe] 


"thin films" 


94.8 


63.4 


50-60 


0.03/0.10 


"thick films" 


92.1 


57.0 


25-30 


0.03/0.11 



strain in the calculations is the same as that found exper- 
imentally. We use a monoclinic angle (3 = 89.5°, consis- 
tent with the experimental data4i All lattice parameters 
are summarized in Table [Q] As reference for the theoret- 
ical lattice parameters we use the values obtained by re- 
laxing the bulk rhombohedral unit cell using U g=2 eV~ 
In the thin films the large in-plane stress leads to a large 
out-of-plane relaxation and a \f2cj a ratio that deviates 
significantly from 1. In the thick films the in-plane stress 
is partially released because of the two different in-plane 
lattice parameters a and b, leading to \/2c/a ratio (or 
\/2c/6 ratio) closer to 1 and nearly no strain in the out- 
of-plane direction. The results for the electric polariza- 
tion and magnetization calculated after relaxing all the 
ionic positions are shown in Table ITTT1 To ensure the in- 
sulating character of the systems we use the LSDA+U 
method for these calculations. 

Due to the larger c/a ratio in the "thin film" struc- 
ture the direction of the polarization rotates further away 
from the [111] direction (towards the [001] direction) 
compared to the "thick film" structure. This leads to 
an increase in the out-of-plane component of the polar- 
ization, in spite of the fact that the effect of strain on the 
absolute value of the electric polarization is rather small. 



Thus, experiments that measure the out-of-plane compo- 
nent of the polarization can expect to see changes in P 
due to this re-orientation effect, even in cases where the 
strain-dependence of the magnitude of the polarization 
is small. 

Comparing the calculated values of the spontaneous 
polarization with the experimental data for (001) ori- 
ented films of different thickness from Ref. |2J (see Ta- 
ble lllip shows that the theoretical values agree reason- 
ably well with the experimental data for the thin films 
(perpendicular component of the polarization). The re- 
ported experimental value for the thicker films seems to 
be smaller than the calculated value. This could be due 
to incomplete switching of the polarization in the thicker 
films. 

In summary, the strong dependency of the polarization 
from the film thickness reported in Ref. is probably 
a sum of two effects: (i) the rotation of the polarization 
away from the [111] direction in the thinner films due to 
the increased c/a ratio, and (ii) incomplete switching of 
the polarization in the thicker films. 



C. Strain dependence of the magnetization 

It can be seen from Tables [|] and IIIII that there is 
no significant strain dependence of the magnetization in 
BiFeOa. The absolute value of the magnetization (re- 
sulting from the canting of the Fe magnetic moments) 
is significantly larger for the traditional J dependent 
LSDA+U treatment compared to the simplified approach 
corresponding to J = eV (see Sec.|nj but in both cases 
this value is not significantly changed by strain, either for 
the rhombohedral symmetry (Table P) or for the mono- 
clinically strained structure (Table |LU|) . 



D. Influence of oxygen vacancies on the electric 
and magnetic properties 

To investigate the influence of oxygen vacancies on the 
ferroelectric and magnetic properties of BiFeOs we per- 
form calculations for a unit cell doubled along one of 
the rhombohedral lattice vectors (so that it contains four 
formula units), in which we remove one oxygen atom 
and then relax all ionic positions in the supercell. We 
then calculate total and partial densities of states, elec- 
tric polarization, and magnetization. We do this for 
the rhombohedral bulk structure as well as for the "thin 
film" monoclinically distorted structure described in Sec- 
tion IIII Bl (see Table [HJ| . The use of a doubled unit cell 
artificially reduces the symmetry of the rhombohedral 
system to the monoclinic space group Bb (in the case 
of the monoclinically distorted structure this is already 
the space group of the original unit cell). In this lowered 
symmetry there are three inequivalent groups of oxygen 
anions. Removing an oxygen anion from one of the three 
groups in turn leads to three different arrangements of 



6 



TABLE IV: Polarization \P a \ and magnetization M s calcu- 
lated within the LSDA+U method for the systems with oxy- 
gen vacancies. The upper (middle) three lines correspond to 
the three vacancy configurations based on the rhombohedral 
bulk structure (monoclinic "thin film" structure). The last 
line corresponds to the tripled unit cell. For the magnetiza- 
tion, the first value refers to U — 2 eV, J = eV and the 
second value refers to U = 3 eV, J = 1 eV (see Sec. [H] for 
details) . 





\P a \ [/iC/cm 2 ] 


Ms \p,B /Fe] 


I 


96.6 


0.01/0.07 


II 


95.0 


0.07/0.14 


III 


97.1 


0.05/0.11 


I 


97.1 


0.02/0.08 


II 


94.6 


0.07/0.14 


III 


99.9 


0.05/0.10 


tripled 


96.5 


0.06/0.12 



oxygen vacancies. In the following we refer to these dif- 
ferent configurations (in a somewhat arbitrary way) as 
configurations I, II, and III. The corresponding supercells 
have a vacancy concentration of 8.3 %, which means that 
one out of 12 oxygen anions is missing. This can formally 
be written as BiFe03_5 with 6 = 0.25. 

All systems are metallic within the LSDA but become 
insulating within the LSDA+U. Table ITVl shows the val- 
ues for the electric polarization and magnetization for all 
supercells calculated within the LSDA+U method. Re- 
markably, the electric polarization is not significantly af- 
fected by the presence of oxygen vacancies and the ac- 
companying structural and electronic distortions. This 
again reflects the high stability of the ferroelectric con- 
figuration in this system. Even the combined presence 
of epitaxial strain and oxygen vacancies in the mono- 
clinically distorted structures does not lead to significant 
changes in the ferroelectric polarization P s . 

Figure[21shows the densities of states (total and partial 
Fe d) for vacancy configuration I based on the rhombo- 
hedral bulk structure together with the corresponding 
densities of states for the same structure without oxygen 
vacancies. The densities of states for the other systems 
all look very similar to that shown in Figure [21 The 
four Fe cations in the supercell containing the vacancy 
can be divided in two classes which we call "Fe 3+ " and 
"Fe 2+ ". The partial d densities of states for the two 
"Fe 3+ " cations are very similar to the case without the 
oxygen vacancies: the (local) majority spin states are 
completely filled and the (local) minority spin states are 
completely empty (apart from a small contribution aris- 
ing from the hybridization with the O 2p states), indi- 
cating a d 5 high-spin electronic configuration. The un- 
occupied minority d states are split into the ti g and e g 
manifolds characteristic of the predominantly octahedral 
symmetry of the crystal field. The partial densities of 
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FIG. 2: Total (gray shaded) and partial Fe d densities of states 
for the unstrained rhombohedral structure (upper panel) and 
for vacancy configuration I (lower panel) calculated for U e s 
= 2 eV. Minority spin states are shown with a negative sign. 
The full lines correspond to the "Fe 3+ " ions, the dotted lines 
in the lower panel corresponds to the "Fe 2+ " ions. The dashed 
vertical lines indicate the highest occupied energy levels. Zero 
energy is set to the upper edge of the "Fe 3+ " majority spin 
bands. 

states for the "Fe 2+ " cations are significantly different 
from this. Here, the t2 g minority states are partially filled 
by approximately one electron (indicating a high spin d 6 
electron configuration) and there is a small gap between 
the occupied and unoccupied minority ti g states. The 
densities of states therefore suggest a picture of distinct 
"Fe 2+ " and "Fe 3+ " cations with d 5 and d 6 electron con- 
figurations respectively. 

Although the presence of both "Fe 2+ " and "Fe 3+ " 
seems rather obvious from the analysis of the densities 
of states, the local charges, obtained by integrating the 
partial densities of states up to the Fermi-energy, differ 
only slightly (~ 0.1 e) for the two types of Fe cations. 
Furthermore, these charges do not represent the for- 
mal charges corresponding to the electron configurations 
mentioned in the previous paragraph. The problem of 
assigning local (static) charges based on the continuous 
electron density in periodic solids is well known (see e.g. 
Ref . I42T) . In practice, local quantities are found by defin- 
ing spheres around the atomic sites and projecting the 
Bloch-functions onto a local basis. Since the radii of such 
spheres are arbitrary within reasonable limits, and the 
intermediate region between these spheres is either ne- 
glected or double-counted, the local quantities defined in 
this way are not uniquely defined. Furthermore, due to 
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FIG. 3: Schematic picture of charge ordering in BiFe02.7sfor 
vacancy configuration I. Large white spheres represent oxy- 
gen vacancies, Large (small) gray spheres represent "Fe 2+ " 
("Fe 3+ ") ions. Different shadings indicate the two different 
antiferromagnetic spin orientations. For simplicity all ions 
are placed at the corresponding ideal positions within the cu- 
bic perovskite structure. On the left side stripes of "Fe 2+ " 
ions along the [110] direction can be identified. On the right 
side the quasi-planar arrangement of the oxygen vacancies is 
visible. 

the formation of band-states within a periodic solid, the 
angular momentum I is no longer a good quantum num- 
ber. Even a pure "d band" , constructed as the Bloch-sum 
of localized d orbitals, usually contains some s aiidp char- 
acter when projected onto a local I basis. It is therefore 
clear that locally defined charges are poor indicators for 
the oxidation states of the ions. The clear qualitative 
differences in the local densities of states appear to be 
much more appropriate for this purpose. 

The three-dimensional arrangement of "Fe 2+ " and 
"Fe 3+ " is shown schematically in Fig. [3] for vacancy con- 
figuration I. The "Fe 2+ " cations appear on the sites ad- 
jacent to the oxygen vacancy. The quasi-planar arrange- 
ment of oxygen vacancies shown in Figure y] is a re- 
sult of the particular restricted supercell geometry. A 
more isotropic distribution of oxygen vacancies and a 
smaller (and probably more realistic) vacancy concen- 
tration would require the use of larger supercells. 

To validate our results also for a slightly smaller va- 
cancy concentration we repeated our calculations for a 
tripled unit cell containing six formula units of BiFe03_a, 
corresponding to a vacancy concentration of 5.6 % or 

5 sa 0.17. In this case we describe the rhombohedral lat- 
tice using hexagonal lattice vectors where one unit cell of 
the hexagonal lattice contains three lattice points of the 
rhombohedral lattice. The resulting hexagonal unit cell 
corresponds to space group P3 and in this case there are 

6 groups of inequivalent oxygen anions. To reduce the 
computational effort, we only treat one possible arrange- 
ment of oxygen vacancies for the tripled unit cell which 
is obtained by arbitrarily removing one of the oxygen an- 
ions. This calculation is done only for the rhombohedral 
bulk lattice parameters and using the LSDA+U method. 

The calculated polarization for the tripled unit cell 
with oxygen vacancy is shown in Table Hvl Also in this 
case the polarization is not significantly affected by the 
presence of the oxygen vacancies. Similar to the doubled 
unit cells, we obtain two distinct classes of Fe cations 



that can be interpreted as Fe 2+ and Fe 3+ , with local Fe 
d densities of states very similar to those shown in Fig- 
ureH In this case the ratio of "Fe 2+ " to "Fe 3+ " is 1:2, as 
required by the charge neutrality of the system. Again, 
the "Fe 2+ " cations appear on the sites adjacent to the 
oxygen vacancy. Therefore, although our supercells en- 
force a rather artificial ordered arrangement of oxygen 
vacancies and the corresponding vacancy concentrations 
are relatively high, we conclude that the incorporation 
of oxygen vacancies in BiFeC>3 leads to the formation of 
Fe 2+ on the sites adjacent to the vacancy. 

We now turn our attention to the effect of vacancies 
on the macroscopic magnetization which is also shown 
in Table IIVI Since the pairs of Fe 2+ cations are always 
situated on neighboring positions of the magnetic lat- 
tice, and are therefore antiferromagnetically aligned, no 
net magnetization results from a ferrimagnetic arrange- 
ment of Fe 2+ and Fe 3+ cations as occurs, for example, in 
magnetite^ 3 " The magnetization is therefore still caused 
entirely by the small canting of the mainly antiferromag- 
netically oriented magnetic moments of the Fe cations, 
as shown in Ref. |2fJ (weak ferromagnetism)m& Two ob- 
servations can be made by inspection of Table llVl First, 
there is no significant difference between the strained and 
unstrained systems. Second, in some cases the magneti- 
zation is enhanced compared to the bulk value, whereas 
in other cases it is decreased. The differences are slightly 
more pronounced for the data set corresponding to U = 
2 eV, J = eV than for the data set corresponding to 
U = 3 eV, J = 1 eV, but no clear trends can be iden- 
tified. The observed changes in magnetization for the 
various vacancy configurations cannot be understood by 
considering simple changes in coordination between the 
magnetically coupled cations (which are the same for all 
vacancy configurations) but rather depend on the details 
of the structural relaxation of both cations and anions. 
This is not surprising since the Dzyaloshinskii-Moriya in- 
teraction which is responsible for the canting of the mag- 
netic moments^ is closely related to the superexchange 
interaction which in turn is known to be very sensitive 
to small structural changes 44 The observed changes in 
magnetization caused by the presence of oxygen vacan- 
cies cannot explain the strong increase in magnetization 
reported for the thin films of BiFeOa^ 



IV. SUMMARY AND CONCLUSIONS 

In summary, our calculations show that the ferroelec- 
tric polarization in multiferroic BiFe03 is extremely in- 
sensitive to both strain and the presence of oxygen va- 
cancies. This insensitivity of the polarization is due to 
a corresponding insensitivity of the ionic displacements 
and is in striking contrast to what has been found for 
most conventional perovskite ferroelectrics. At present it 
is not clear if this stability is due to a general high stabil- 
ity of the ferroelectric state in BiFeC>3, reflected also in 
its large ionic displacements and high Curie temperature 
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(T<7 = 1123 K), 37 if this stability is a special feature of the 
different mechanism driving the ferroelectric distortion in 
this class of Bi-containing multiferroic materials, namely 
the stereochemically active Bi lone electron pairj2& or if 
it is due to the special geometry of the i?3c space group 
containing oxygen octahedra rotations in addition to the 
polar displacements. Future studies of the strain depen- 
dence of the ferroelectric polarization in the high Tq fer- 
roelectric LiNb03 and the lone-pair active multiferroic 
BiMnC>3 will help to solve this issue. 

The incorporation of oxygen vacancies in BiFeC>3 leads 
to the formation of Fe 2+ which can be identified by 
the clear qualitative differences in the local densities of 
states, but the actual charge disproportionation is small. 
The presence of oxygen vacancies can affect the value of 
the macroscopic magnetization, although the observed 
changes are too small to explain the strong increase of 
the magnetization reported for thin BiFeC>3 film, 23 These 
effects of oxygen vacancies are independent of the strain 
state of the system. 

Recently, measurements of the electric polarization for 
(001), (101), and (111) oriented films of BiFe0 3 (200 nm 
thickness) were reported 2 ^ and the experimental data 
could be explained as resulting from different projections 
of the same polarization vector. This supports the notion 
that both magnitude and orientation of the polarization 
in BiFeC>3 are not significantly affected by strain, since 
the strain tensor differs considerably for different film ori- 
entations. The reported value of ~ 100 /iC/cm 2 for the 
(111) oriented filmSi agrees well with the polarization 
values calculated in this work. 



Here, n^ m , are the elements of the orbital density matrix 
and n — ^2 sm nf£ m is the total number of d electrons at 
the corresponding ion. Summation over all sites contain- 
ing d electrons has been suppressed for simplicity. 

The corresponding expression for the LSDA+U ap- 
proach of Anisimov, Liechtenstein and coworkers is>2i 



E — i?LSDA 
1 



S2 

37714 



{m,s} 

- (mim3|Fee|m 4 m3) n^n^} 



2 -„(„-l) + ^2J »>*-!) • (A2) 



Here, n s = J2 m n mm i s the total number of d elec- 
trons with spin s and {m\m 3 \V ee \rn 2 m,4) are the matrix 
elements of the screened electron-electron interaction. 
These matrix elements are defined in terms of Slater- 
integrals F k as: 

{mim z \V cc \m 2 m i ) = 2J a* (mi, m 3 , m 2 , m 4 ) F k , 

k 

(A3) 

with 



a k (m 1 ,m 3 ,m 2 ,m i ) 



4tt 



2k- 



- J2 ( lm i\ Y k q \lm 2 )(lm 3 \Y k * q \l mi ) . (A4) 
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APPENDIX A: EQUIVALENCE OF THE TWO 
LSDA+U APPROACHES FOR J = 

In the LSDA+U approach of Dudarev et alM* the total 
energy of the system is expressed as follows (Equation (5) 
of Ref. I32L generalized for noncollinear spin systems): 



E = -Elsda + ~y I " 



E 



T! SS ,77 S f 

mm mm 



(Al) 



The Slater integrals are related to the parameters U and 
J via U = F°, J — (F 2 + F 4 )/U, where one usually 
sets F 2 /F 4: = 0.625. For J = the only nonvanishing 
term in the sum of Eq. (|A3|) is for k = 0, and one can see 
from Eq. (|A4|) that a (mi, m 2 , m 3 , m 4 ) = 5 mi „ l2 8 m3mi . 
For J = the matrix elements (|A3|) are therefore simply 
given by: 



(mim 3 \V ee \m 2 m4) 
Using HA5f) in (|A2|) one obtains: 



E — -Elsda 




E 



71 SS ,77 S f 
mm mm 



(A5) 



(A6) 



m.m's.s' 



which is just Eq. (XT} with U cfi = U. The LSDA+U 
approach of Dudarev et alS^ is therefore included in the 
approach of Ref. |3l] as the special case J = 0. 



m.m' ,s,s' 
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